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The theory of heat transfer by electromagnetic radiation is based on the radiative transfer equation 
(RTE) for the radiation intensity, or equivalently on the Boltzmann transport equation (BTE) for the 
photon distribution. We focus in this review article, after a brief overview on different solution methods, 
on a recently introduced approach based on truncated moment expansion. Due to the linearity of the 
underlying BTE, the appropriate closure of the system of moment equations is entropy production 
rate minimization. This closure provides a distribution function and the associated effective transport 
coefficients, like mean absorption coefficients and the Eddington factor, for a general number of moments. 
The moment approach is finally illustrated with an application of the two-moment equations to an 
electrical arc. 



PACS: 44.40.+a, 52.25.Os, 95.30.Jx 



1. Introduction 

Heat radiation refers to electromagnetic radiation emitted by thermally excited degrees of 
freedom of matter. If both the matter and the radiation field are at thermodynamic equi- 
librium, well-known relations from thermodynamics exist between the temperature T, the 
characteristic radiation frequency v, the energy density and the pressure p rac i of the 

radiation field. These are Wien's displacement law, v = 5.88 ■ 10 10 T (v in Hz and T in K), 
the caloric equation of state, E") = 7.57 ■ 10~ 16 T 4 (in J/m 3 ), and the thermal (or thermo- 
dynamic) equation of state, p rac ( = E^) /3. It is then straight-forward to derive the Stefan- 
Boltzmann law for the po wer emitted by a black body, Q = 5.67 ■ 1CP 8 T 4 in units of W/ m 2 (cf. 
ILandau & Lifshitzl ^20051) 1. In typical applications, heat radiation is relevant in the frequency 
range of 10 11 — 10 16 Hz, including the upper part of the microwave band, the infrared, the 
visible light, and the lower part of the ultra-violet band. 

In many cases, be it for engineering purposes like electric arc radiation modelling, or related 
fundamental scientific problems like in stellar physics, radiation is usually not at thermal 
equilibrium. The present chapter of this book aims to give a focused overview on the the- 
ory of radiative heat transfer, i.e., energy transport by heat radiation that can be in a gen- 



eral nonequilibrium state, while matter is at local thermodynamic equilibrium. With em- 
phasis on models based on partial differential equations for the radiation energy density, 
heat flux, and (if necessary) higher order moments, we will particularly discuss a powerful 
method for the determination of effecti ve transport coefficients, which has been recently de- 
ve loped bvlChristen & Kassubekl j2009l). General mon ographs on rad iative transfer are given 
by lChandrasekharl jl960|) . ISiegel & Howelll jl992h . and lModestl i2003h . to mention a few. 
In Sect. |2 the basic definitions and equations for radiative heat transfer will be introduced. 
There are two equivalent descriptions of radiation, either in terms of the specific radiation in- 
tensity (or radiance), I v (x,Cl,t), or the photon distribution function, f v (x,Cl,t). Here, t, x, v, and 
O denote time, position, frequency, and direction (normalized wave-vector), respectively. Fre- 
quency dependence will always be indicated by an index v. The associated transport equa- 
tions for I v and f v are named the radiative transfer equation (RTE) and the Boltzmann transport 
equation (BTE), respectively. The number of photons in the volume element d 3 x at position x 
and time t, in the frequency band dv at v, and in direction CI within the solid angle dO, equals 
f v {x,Cl,t)d i xdvdCl. The intensity is then given by l v (x,Cl,t) = chvf v , where hv is the photon 
energy, h = 6.626 ■ 10~ 34 [s is Planck's constant, and c = 2.998 ■ 10 8 m/s is the vacuum velocity 
of light (cf.He3ll968j)). I v is the energy current density per solid angle in direction O. 
The RTE (or BTE) is an integro-differential equation for I v (or f v ) in the 6-dimensional phase 
space corresponding to position, frequency, and direction, and describes the temporal change 
of l v (or /,,) due to emission, absorption, and scattering by the matter. Finding an appropri- 
ate solution is generally a highly sophisticated task, and can be significantly impeded by a 
complicated frequency dependence of the radiation-matter interaction. Moreover, radiation 
problems in science and engineering often require a self-consistent solution of the coupled 
equations for radiation and matter. For instance, a treatment of radiation in hot gases or 
plasma involves, besides the RTE, the gas-dynamic balance equations for mass, momentum, 
and energy (or temperature). Despite of the recent year huge progress in computational tech- 
nologies, an exact solution of the complete set of coupled equations is still unfeasible, except 
for some especially simple cases. As a consequence, in the course of time a number of meth- 
ods for approximate solutions of the RTE have been developed. In Sect. [3j we will therefore 
briefly discuss a selected list of important approximation concepts. Methods based on trun- 
cated momentum expansions will be emphasized, and the need of a reliable closure method 
for the determination of the transport coefficients occurring in these equations will be moti- 
vated. 

In Sect. |4]we will argue that a recently introduced approach for the closure based on entropy 
production rate is superior to other closures us ed up to date . The theory of radiation in ther- 
mal equilibrium dates back to seminal work by lPlanckl (|l906h . In chapter 5 of his book Planck 
emphasizes that photons, unlike a normal gas of massive particles, do not interact among 
themselves but interaction with matter is needed for a relaxation to the thermal equilibrium 
state. As is often the case in many applications of radiative transport, we will assume that the 
medium, be it condensed matter, gas or plasma, is at local thermodynamic equilibrium (LTE) 
and can thus be described locally by thermodynamic quantities like temperature, chemical 
potential, and the same. It is then the equilibration process of the photon gas to this LTE state 
that determines the details of the heat transfer by radiation in the medium. As is well-known 
from thermodynamics, equilibration is related to entropy production, w hich plays an im por- 
tant role in understanding the behavior of nonequilibrium radiation (cf . lOxeniu si <1966fr and 
iKrollI dl967l) '). In fact, various authors have shown that the state of radiation is often related to 
optima of the entropy production rate. Whether the optimum is a maximum or a minimum, 
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depends on the specific details of the system under consideration, particularly on convexity 
prope rties of the optimization problem, and particularly, the constraints. For instance, lEssexl 
jl984h has shown that the entr opy producti on rate is minimal in a gray atmosphere in local 
radiative equilibr i um. L ater on lEssexl fi 997) applied his approach also to neutrino radiation. 
IWurfel & Ruppell Jl985h and lKabelad J1994I) discussed entropy production rate maximization 
by intr oducing an effective c hemical potential of the photons, related to their interaction with 
matter. ISantillan" etal.ljl998h showed that for a constraint of fixed radiation power, black bod- 
ies maximize the entropy production rate. 

The underlying r eason for the success of entropy production rate principles has been recog- 
nized already by lKohle"rl fl948h , who has shown that the stationary solution of the BTE that 
is linearized at the equilibrium distribution, satisfies a variational principle for the entropy 
productio n rate. Kohler 's principle has b een widely used to determine linear transport coeffi- 
cients (cf. IZimanl il956l) and refs. cited in lMartyushevI j2006l) ). The important property of the 
RTE (or the BTE for photons) is its linearity over the whole nonequilibrium range, provided 
the interaction with the LTE-medium consists of single-photon processes only. This linearity is 
thus not an approximation as it was in Kohler 's work, but holds for arbitrarily large deviations 
from thermal equilibrium of the photon gas. The absence of interaction between photons is 
thus the reason for the success of the concept beyond small deviations of f v from equilibrium. 
Consequently, the entropy production rate is the appropriate basis for the determination of 
the nonequilibrium distribution l v (or /„) and the effective transport coefficients for radiative 
heat transfer in the framework of a truncated moment expansion. In Sect. [5] the transports 
coefficients, i.e., the effective absorption constants and the Eddington factor, are calculated 
for some specific examples. A practical reason for selecting moment equations for modelling 
radiative transfer is the convenience of having a set of structurally similar equations for the 
simulation of the complete radiation-hydrodynamics problem. Both the hydrodynamic equa- 
tions for matter and the moment equations for radiation are hyperbolic partial differential 
equations and can thus be solved on the same footing. In Sect. I5.4l gives some remarks on the 
requirement of hyperbolicity. For numerical simulations boundary conditions must be spec- 
ified, these will be discussed in Sect. [6] Finally, Sect. [7] will then provide some simulation 
results for a simplified example of electric arc radiation. 

2. Basics of Radiative Heat Transfer in Matter 

The radiation intensity, I v (CI), is governed by the radiative transfer equation (RTE), 

^d t i v + ci-vi v = K V (B v -i v ) + a v ^J s2 dnp v (n,ti)i v {Ci)-h}j , (i) 

which has to be solved in a spatial region defined by the physical problem under consider- 
ation. Phase coherence and interference effects are disregarded when considering thermal 
radiation, and we will also not consider polarization effects. The BTE is simply obtained by a 
replacement of I v by f v in the RTE. The left hand side gives the total rate of change of I v (Cl), 
divided by c, along the propagation direction O. This change must be equal to the expression 
on the right hand side, which consists of a sum of specific source and sink terms due to the 
radiation-matter interaction. In the absence of any interaction, e.g., in vacuum, the right hand 
side vanishes, which describes the so-called (free) streaming limit associated with a radiation 
beam, or the ballistic propagation of the photons. In the presence of interaction, however, pho- 
tons are generated by emission and annihilated by absorption, described by k v B v and —k v I v , 
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respectively. Here, B v is the Planck function for thermal equilibrium, 

B v = ¥»W , (2) 

where 

J e l) - I (3) 

exp(hv/k B T)-l ^ 

is the Bose-Einstein distribution for thermal equilibrium photons (cf. lLandau & Lifshitzl 
j2005t) ) with the Boltzmann constant k B = 1.381 ■ 10~ 23 J/K and the local temperature T = T(x) 
of the LTE medium. The coefficient k v is the macroscopic spectral absorption coefficient 
in units of 1/m, and is generally a sum of products of particle densities, absorption cross- 
sections, factors [1 — exp(— hv/kgT)], and depends thus not only on frequency but also on the 
partial pressures of the present species and the temperature. Often, opacities referring to K v /p 
are discussed in the literature, where p is the mass den sity of the m atter. The macroscopic k v 
includes spontaneous as well as induced emission (cf. iTienl ll968h ). Additionally to inelastic 
absorption-emission processes, Eq. 10 includes elastic (or so-called coherent or conservative) 
scattering. Incoming photons of frequency v from all directions fi are scattered with probabil- 
ity p v (n,n) into direction fi. It is assumed that the phase-function p v (Ci, Ci) obeys symmetry 
relations a ssociated with recipro city, depends only on the cosine between the directions fi 
and n (cf. IChandrasekharl (1960|V ). and is normalized, (An) 1 J s2 dClp v (Cl,Cl) = 1. Here, the 
Q-Integration extends over S 2 , which denotes the unit sphere associated with the full solid 
angle. The strength of the scattering process is quantified by the spectral scattering coefficient 
o~ v in units of 1/m. The ratio o~ v / (k v + a v ) gives the probability that a collision event is a scat- 
tering process, and is sometimes called the (single-scattering) albedo. The mean free path of the 
photons is the inverse of k v + a v . 

Because the bracket proportional to o~ v in Eq. {TJ vanishes for fi-independent l v , the RTE can 
be written in the simple form 

-dtlv + Cl- VI V = £(B V -I V ) , (4) 

where the linear, self-adjoint, positive semi-definit^] operator C is defined by the right hand 
side of Eq. {T} and consists of an algebraic term and an integral term. 

The RTE has to be solved with appropriate initial conditions, I v (x,Cl,t = 0), and boundary 
conditions on the surface of the spatial domain under consideration. Because the RTE is a first 
order differential equation, the determination of each ray requires the knowledge of I v (O) on 
the domain surface and in directions O pointing into the domain. The behavior of the bound- 
ary is characterized by the radiation it emits, and the way it reflects impinging radiation. If 
one denotes the emittance of the boundary at position x w by e(x w ), the reflectivity by r(x w ), 
and the normal vector of the boundary surface by n(x w ), the boundary condition generally 
reads (cf . iModestl j20O$l ) 

I v (\ w ,Cl,t) = e(x w )B,/(x w ) + / dtl I n(x w ) ■ O I r(x w ,Cl,fl)I v (x w ,Cl,t) . (5) 

7n(x w )n<o 

The integration runs over all n associated with radiation coming from the bulk domain 
towards the surface, while CI is pointing into the domain. For a smooth surface where a 



1 Note that a negative eigenvalue would immediately lead to an instability. 
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normal vector n(x w ) can be defined, this solid angle corresponds to half of the sphere S 2 . 
This general boundary condition can be simplified for special limit cases. For instance, a 
black surface has r = and e = 1, a diffusively reflecting surface has r(x w ,fi,f)) = r(x w ) / n, 
and a specularly reflecting surface has r(x w ,n,0) oc 5{Cl s — CI), where Cl s = CI — 2(fi • n)n 
is the direction from which the ray must hit the surface in order to travel into the direction of 
O after specular reflection. 

We conclude this section by listing the basic equations for the LTE matter to which radiation 
is coupled. In general, LTE implies that at each point in space, the caloric and thermody- 
namic equations of state are locally valid. The respective equations relate the specific energy 
e = e(p,T) and the pressure p = p{p,T) to the mass density p and the temperature T of the 
matter. The spatio-temporal dynamics of the thermodynamic variables and, if relevant, the 
flow velocity u, is then given by the hydrodynamic balance equations for mass, momentum, 
and energy. For a single component (non-relativistic) medium 



where n ma t, j e , and etot = £ + u 2 /2 are the momentum stress tensor, the energy flow density, 
and the total energy density. Together with the equations of state, Eqs. |6)-|8) constitute seven 
equations for the seven variables p, p, T, e, and u. The right hand sides, p, f, and W are the 
mass source density, the force density, and the heat power density, respectively. The effect 
of radiation on matter may occur in these three source terms. Fo r instance, a m ass source 
may appear at a solid wall due to ablatio n by radiation (see, e.g. | Christenl J2007|) ). and the 
radiation pressure may act as a force (cf . iMihalas & Miha las ( 1984)). These two effects are 
often negligible in e ngineering applications or are i mportant only in special cases, like ablation 
arcs as discussed by lNordborg & Iordanidisl (|2008h . However, the heat exchange described by 
W can in general not be disregarded, and will play an important role in the theory below. 
The back-coupling of the matter on radiation, as mentioned before, occurs in the expressions 
on the right hand side of Eq. {T}, which depend gene rally on p (or p) and T. A n extensive 
monograph on radiation hydrodynamics is provided by lMihalas & Mihalasl dl984l) , and a short 
introd uction that fits well to the present chapter is given by the lecture notes of IPomraningl 
<1982h . 

3. Approximation Methods 

The extreme difficulties to solve the RTE exactly for real systems caused the development of 
various approximation methods. There are two additional reasons for the use of approxi- 
mations. First, in many cases the behavior of the matter is of interest, while it is sufficient to 
consider the radiation as a means of (nonlocal) interaction; hence only the radiative heat flux is 
needed, which enters the power balance equation for the matter via the heat power density W. 
As W equals the negative divergence of the radiation energy flux density, a radiation model 
would be convenient that is confined to this flux and to the lower order moments, which is 
here a single one, namely the radiation energy density. Secondly, radiation often behaves in 
two different specific ways. In a transparent medium absorption and scattering are weak, and 
radiation propagates as beams; full absence of interaction with matter refers to the so-called 
free streaming limit. In an opaque medium, on the other hand, absorption, emission and /or 
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scattering is strong, and the radiation diffuses isotropically. In the extreme diffusive limit the 
Rosseland diffusion approximation applies, where radiative transfer is modelled by an effective 

heat conductivity of the matter (cf . ISieg el & Howelll Jl992h ) given by 16a SB T 3 /3a ( / i] . Here 

a SB = 2n 5 k 4 B /15h 3 c 2 = 5.67 ■ 1CT 8 W/m 2 K 4 is the Stefan-Boltzmann constant and a^ S) is the 
Rosseland mean absorption to be discussed later. For the two behaviors a ballistic (beam-like) 
and a diffusive description, respectively, are the appropriate 'zero order' models with effec- 
tive transport coefficients, and deviations from the limits may be treated by small corrections. 
Models based on one of these two limit cases can strongly reduce the computational effort. 
However, in many real systems radiative transfer is in between these limits such that more 
sophisticated methods must be involved. 

In the following, a short list of some relevant approximation methods is given. The selection 
is n ot complete , as other approaches exist, like ray tracing and radiosity-irradiosity meth- 
ods teevl(l2006l)'). or so me rather heuristic method s like the Pj/3-approximation discussed by 
lOlson et al I bOOCh andU immons & Mihalasl J2000h . Furthermore, we will not discuss the is- 
sue of discretization methods concerning position space like finite differences, volume s, or 
elements; although this field would require special recognition (cf. lArridge et alJ teOOOl) and 
Refs. cited therein) it is beyond the purpose of this chapter. Needless to say that there is not 
a unique best method but every approach has its advantages and disadvantages for practical 
use, and the appropriate choice depe nds usually on the problem under consideration. Exh aus- 
tive overviews can be found, e.g., in lDuderstadt & Martini Jl979t l . [Siegel & Howelll 1 19921) and 
literature cited in the following three subsections. Subsequently, we will then focus in subsec- 
tion [3T4] on approximations based on moment expansions, and particularly on the closure of 
the moment equations that will be discussed in Sect. [4] 

3.1 Net Emission 

The net emission approximation is probably the most simplistic radiation model. It assumes 
a semi-empirical function W(T,p,£) in Eq. (8). Additionally to temperature and pressure, 
it depends on parameters £ of the radiating object. It is somet imes used, fo r instance, in 
comp u tational fluid dyna mics simulations of electrical arcs (cf . iLowkd jl970h , IZhang et all 
Jl987h . ISeeger et alJ J2006l) ), where the only parameter £ is the arc radius. Although such a 
description is very convenient in numerical simulations and sometimes even provides useful 
results, it is obviously oversimplifying and without any rigor. Furthermore, reliable accuracy 
requires, for the determination of the function W(T,p,£), a parameter study based on a more 
fundamental radiation model or on elaborate experiments. 

3.2 Monte Carlo 

Mo nte Carlo simulations refer to r andom sampling methods (see, for instance lYang et all il995l) 
and iDuderstadt & Martini il979l) ), which are based on computer simulations of a number of 
photons. Their deterministic dynamics corresponds to the ballistic motion with speed of light. 
Emission, absorption, and scattering processes are simulated in a probabilistic way by appro- 
priately determined random numbers for the various processes. Those include, of course, the 
interaction with boundaries of the spatial domain. Final results, like the radiation intensity, 
are determined by averages over many particles. The Monte Carlo co ncept is rather sim ple, 
which leads to a number of advantages of this method, as discussed by lYang etalHl995l) . Ef- 
ficient applications make u se of specifically imp rove d schemes like imp licit Monte Carlo or 
special versions thereof (cf. iBrooks & FleckMl986h and lBrooks et all J2005t) ). 
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3.3 Discrete Ordinates 

The discrete ordinates method (DOM) considers a finite number of rays passing at every (dis- 
crete) space point. If a number Njj of direction vectors fi; c , k = 1, ...,Nrj> is selected, one has 

l v = ly' S(Ci — f\), such that a set of Nu partially coupled RTE-equations for the differ- 
ent directions and frequencies must be solved. The right hand side of these equations, say 

C(B V — ly ), contains not an integral as Eq. Jl) but a weighted sum. As reasonable minimum 
values for Nq in 3-dimensional realistic geometries are of the order of 10, the computational 
effort is still large. For too small N-q an artifact called "ray effect" may occur, referring to 
spatial oscillations in the energy density. Another error known as "false scattering" or "false 
diffusion", is due to the discretizati on of position space and is linked in a certain way to the 
ray effect as discussed in lReyl 120061) . 

Some further developments based on DOM exist, which make use of a decomposition and dis- 
cretization of the angular space into a finite set of directions, i.e. a finite p artition of the unit 
sphere 5 2 . The methods of partial characteristics iAubrecht & Lowkel (T99J)) and of partial mo- 
ments {Frank et al.l te006|) 1 are examples, the latter being mentioned again in the next section. 
Last but not least, we mention that is h as been proven that the DO M is equivalent, un der cer- 
tain conditions, to the P-N method (cf. iBarichello & Siewertl jl998h and lCullenl l200ll» , which 
is a special kind of the moment approximations to be discussed in the next subsection. 

3.4 Moment Expansions 

Radiation modelling in terms of moments of the distribution l v (or f v ) is convenient because 
the radiation is coupled to the LTE matter in Eqs. {6j-{SJ via the first thr ee (angular) mom ents. 
Moment expansi ons can be formulated in a rather general manner (cf. iLevermorg 1 19961) and 
IStruchtrupl jl998h ). In the following, we define moments based on I v bj0 




with (fi : fi)jy = f\0/. The last line indicates that an infinite number of moments exist in 
general. E„, F v , and Tl v are, respectively, the monochromatic energy density, radiative flux, 
and stress or pressure tensor of the radiation. For convenience, the pref actor (c _1 ) is chosen 
in all definitions such that the moments have the same units of a spectral energy density. 
Similarly, E, F, and n are the spectrally integrated energy density, radiative flux, and pressure 
tensor. In the present units F has the meaning of energy density associated with the average 
directed motion of the photons, and E of the total energy density composed of directed and 
thermal fluctuation parts. Hence, F =| F | < E, which will be important below. 
In thermal equilibrium all fluxes vanish. Then F^ = 0, the stress tensor is proportional to the 
unit tensor with diagonal elements EW/3, and the energy density is given by 

E (eq) = r AnAvBv = ^SB T 4 . (12) 
Jo c 

2 We mention that a moment corresponding to the photon number (obtained by integration over f v ) 
does not appear, partly because the photon number is not a conserved quantity. 



7 



The purpose of a moment expansion is to derive from the RTE or BTE balance equations for 
the moments, either for each frequency v, or for groups of frequencies or frequency bands, or 
for the full, integrated spectral range. Multiplication of the RTE with products and / or powers 
of fijt's, and integration over the solid angle gives for the moments E v , F v , etc. 



etc., where only the first two equations are listed for convenience, but the list still contains 
an infinite number for all moments and for all frequencies. Practical usability calls then for 
a two-fold approximation. First, the list of moments, and thus moment equations should be 
truncated by considering only the N first moment equations. Secondly, the frequency space 
should be discretized or partitioned in some way, in order to end up with a finite set. If the 
spectrum allows a division into a number of well defined frequency bands with approximately 
constant k v and cr v , or a grouping of different frequencies together according to similar values 
of k v and a v , one can average the equations over such partitions. The associated methods 
are sometime s named multi- g roup, multi-band, or multi-bin methods. For details , we refer 
the reader to iTurpaultl <2005h . iRipoll & Wra"y1 <2008h . iNordborg & Iordanidid l|2008h . and the 
literature cited therein. In the following we will consider the equations for the spectrally 
averaged quantities, which are obtained by integration of Eqs. {13}, 114t . etc., over frequency 



etc., where the right hand sides define Pe and Pp, etc. These quantities are still functionals of 
the unknown function I v . All moments, on the other hand, are variables that are determined 
by the full (still infinite) set of partial differential equations, provided reasonable initial and 
boundary conditions are given. 

Now we perform a truncation by using only the first N moment equations. The first N mo- 
ments would then be determined by the solution of these equations, if the right hand sides 
(Pe, Pf/ etc) and the N + l'th moment were known. In the following section we will dis- 
cuss closure methods that determine these unknowns that are supposed to be functions of the 
N moments. Prior, however, we remark that instead of using products of Cartesian coordi- 
nates of fi, one may equivalently consider a representation in terms of spherical coordinates 
(0,cp). The radiation density is then expanded in spherical har monics Y™ (6, tp) ■ If trun cated, 
this approximation corresponds t o the P-N approximatio n (cf. Siegel & Howelll |l992)). The 
prominent P-l approximation (cf. Siegel & Howelll il992(l '). for instance, refers to a truncation 
of the Eqs. Jl3t and fl4t (or Eqs. il5t and fl6t ) after the second equation and considers an 

isotropic Yl v (or IT) with diagonal elements equal to E v /3 (or E /3). 

We also mention again the partial moment approximation (cf. iFrank et"ai1 teOOtj) ). where the 
approaches of DOM and moment expansion are combined in a smart way. As the DOM dis- 
cretizes the angular space in different directions, the partial moment method selects partitions 

A of the unit sphere S 2 and defines partial moments E v "^\ F v ~^\ H v ~^\ etc, where the solid 
angle integration is performed only over A instead of the whole S 2 . The most simple but 
nontrivial partial moment model refers to the forward and backward traveling waves in a 




(13) 



(14) 




(15) 



(16) 
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one-dimensional position space, where the integration occurs over the two half-spheres as- 
sociated with forward and backward directions. According to iFrankl J2007h , this method is 
able to resolve a shock-wave artifact occurring for counter-propagating and interpenetrating 
radiation beams. 



4. Closure Approaches 

The quality of the moment approximation depends on the number of moments taken into 
account, and on the specific closure concept. A closure of a truncated moment expansion re- 
quires in principle knowledge of I v . A simplification occurs if k v , a v , and p v are assumed to 
be constant (gray matter). The right hand sides of the moment equations strongly simplify as 
they can be directly expressed in terms of these constants and linear expressions of the mo- 
ments. But in general matter is non-gray and the absorption and scattering spectra can be 
extremely complex. Furthermore, the N + l'th moment remains still unknown even for gray 
matter. In the sequel we will discuss a few practically relevant closure methods. We will then 
argue that the preferred closure is given by an entropy production principle. 
For clarity we will consider the two-moment example; generalization to an arbitrary number 
of moments is straight-forward. The appropriate number of moments is influenced by the 
geometry and the optical density of the matter. For symmetric geometries, like plane, cylin- 
drical, or spherical symmetry less moments are needed than for complex arrangements with 
shadowing corners, slits and the same. For optically dense matter, the photons behave diffu- 
sive, which can be modelled well by a low number of moments, as will be discussed below. 
For transparent media, beams, or even several beams that might cross and interpenetrate, may 
occur, which makes higher order or multipole moments necessary. 

4.1 Two-Moment Example 

The unknowns are Pj, Pp, and n, which may be functions of the two moments E and F. For 
convenience, we will write 

P E = 4 eff) (E^ - E) , (17) 
P F = -4 6ff) F . (18) 

where we introduced the effective absorption coefficients Kt; and Kp) that are generally 
functions of E and F. Because the second rank tensor n depends only on the scalar E and the 
vector F, by symmetry reason it can be written in the form 

n„,„ = E^ + ^^) , (19) 

where the variable Eddington factor (VEF) x is a function of E and F. Assuming that the under- 
lying space is isotropic, Kj f , and % can be expressed as functions of E and 

v=\, (20) 

with F = F | . Obviously it holds < v < 1, with v = 1 correspondi ng to a fully directed 
radiation beam (free streaming limit). According to IPomranin gl |l982h . the additional E de- 
pendence of suggested or derived VEFs often appears via an effective E-dependent single 
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scattering albedo, which equals, e.g. for gray matter, (kE^ + o~E) / (k + o~)E. 
The task of a closure is thus to determine the effective transport coefficients, i.e., effective mean ab- 
sorption coefficients k^ S \ k^ S \ and the VEF x of E and F (or v). This task is of high relevance in 
various scientific fields, from terrestrial atmosphere physics and astrophysics to engineering 
plasma physics. 

4.2 Exact Limits and Interpolations 

In limit cases of strongly opaque and strongly transparent matter, analytical expressions for 
the effective absorption coeffici ents are often used, which can be determined in princ iple from 
basic gas properties (see, e.g., lAbuRomia & Tier] jl967|) and iFuss & Haminsl (20021^ . In an 
optically dense medium radiation behaves diffusive and isotropic, and is near equilibrium 
with respect to LTE-matter. The effec tive absorption coeffici ents are given by the so-called 
Rosseland average or Rosseland mean (cf . ISiegel & Howdl \\99^\ ) 

K {eii) = ( Kv ) Ro := !o dv ^ n( ' q) t (21) 
Jg 00 dvv^Ky X 'b v rVy' ] ' 

where d v denotes differentiation with respect to frequency, and 

4 6ff) = <*v + <r v )so ■ (22) 

The Rosseland mean is an average of inverse rates, i.e., of times, and must thus be associated 
with consecutive processes. A hand-waving explanation is based on the strong mixing be- 
tween different frequency modes by the many absorption-emission processes in the optically 
dense medium, due to the short photon mean free path. 

Isotropy of II implies for the Eddington factor % = 1/3. Indeed, because YJ^-kk = E> one nas 
then 0^; = SfciE/3, where 5y (= if k ^ I and 5^ = 1 if k = I) is the Kronecker delta. With these 
stipulations, Eqs. (15} and (16) are completely defined and can be solved. 
In a strongly scattering medium (cr v 2> k v ), where F relaxes quickly to its quasi-steady state, 

one may further assume F = — V£/3Kp for appropriate time scales. Hence Eq. (TBI becomes 



1„ „ „ / VE 

c 



f^S) )=4 eff) (£ ( 



which has the form of a reaction-diffusion equation. For engineering applications, E often 
relaxes much faster than all other hydrodynamic modes of the matter, such that the time 
derivative of Eq. (23) can be disregarded by assuming full quasi-steady state of the radi- 
ation. Equation ((23) is then equivalent to an effective steady state gray-gas P-l approximation. 

For transparent media, in which the radiation beam interacts weakly with the matter, the 
Planck average is often used, 

J-dvv^ 

In contrast to the Rosseland mean, the Planck mean averages the rates and can thus be 
associated with parallel processes, because scattering is weak and there is low mixing 
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between different frequency modes. In contrast to the Rosseland average, the Planck average 
is dominated by the largest values of the rates. Although in this case radiation is generally not 
isotropic, there are special cases where an isotropic FT can be justified; an example discussed 
below is the v — > limit in the emission limit E /EW — > 0. But note that x = 1 often occurs in 
transparent media, and consideration of the VEF is necessary. 

In the general case of intermediate situations between opaque and transparent media, heuris- 
tic interpolations between fully diffusive and beam radiation are somet imes perform ed. Ef- 
fective absorpt ion coefficients have been constructed heuristically by iPatchl fl967h . or by 
ISampsonl jl965h by interpolating Rosseland and Planck averages. 

The consideration of the correct stress tensor is even more relevant, because the simple x — 1 /3 
assumption can lead to the physical inconsistency v > 1. A common method to solve this 
problem is the introduction of flux Iimiters in diffusion a pproximations, where the effec - 
tive diffusion con stant i s assumed to be st ate-dependent (cf . iLevermore & Pomrarungl Jl98lh , 
iPomraningl jl98ll) . and iLevermord jl984) , and Refs. cited therein). A similar approach in 
the two-moment model is the use of a heuristically constructed VEF. A simple class of flux- 
limiting VEFs is given by 

i + zvi 

X = —, (25) 

with positive ;'. These VEFs depen d only on v , but not additionall y separately on E. The 
cases 7 = 1 and j = 2 are attributed to lAuerl Jl984l) and lKershawl dl976h . respectively. While the 
former strongly simplifies the moment equations by making them piecewise linear, the latter 
fits quite well to realistic Eddington factors, particularly for gray matter. 

4.3 Maximum Entropy Closure 

An often used closure is based o n entr opy maximization (c f . jMinerbol l |l978h . lAnile et all jl99lh . 
ICernohorsky & Bludmanl dl994f) , and iRipoll et alj teOOlh 'lFl This closure considers the local 
radiat ion entropy as a function a l of l v . The entr opy i s defined at each position x and is given 
by (cf. lLandau & Lifshid fcOOSh . lOxeniusI jl966h . and|Krol]| <|l967h ^ 

f 2v 2 

S r ad[A'] = —kg / dCldv—^- (n V /lnn v , — (1 + w v )ln(l + n v )), (26) 

where 

Mx ' n) =S^ (27) 

is the photon distribution for the state (v,fi)@ At equilibrium [27} is given by @. l v is then 
determined by maximizing S rat j [I v ], subject to the constraints of fixed moments given by Eqs. 
{9), ( TTOt etc. This provides l v as a function of v, Q, E and F. If restricted to the two-moment 
approximation, the approach is so metimes called the M-l closure. It is generally applicable t o 
multigroup or multi band models jCullen & Pomranind jl98d)jRipolil <2004 . lTurpaultl j2005h , 
IRipoll & Wrayl <2008h ) and partial moments <Frank et al1<2006l).lFrank| j2007h ^ as well as for an 
arbitrarily large number of (generalized) moments jStruchtrupl <1998f)1. It is clear that this clo- 
sure can equally be applied to particles obeying Fermi statistics (see lCernohorsky & Bludmanl 

3 In part of the more mathematically oriented literature, the entropy is defined with different sign and 
the principle is called "minimum entropy closure". 

4 Note the simplified notation of a single integral symbol j in Eq. {26) and in the following, which is to 
be associated with full frequency and angular space. 
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<1994h and lAnileetallteOOOk 

Advantages of the maximum entropy closure a re the mathematic al si mplicity and the miti- 
gation of fundamental physical inconsistencies dLevermord il996|) and|Frank (2007j)). In par- 
ticular, there is a natural flux limitation by yielding a VEF with correct limit behavior in both 
isotropic radiation (x — > 1/3) and free streaming limit (x — > 1): 



XM^ = l~l^--^ (28) 

that depends only on v. Furthermore, because the o ptimization prob lem is conve^fl, the 
uniqueness of the solution is ensured and, as shown by lLevermorel (|l996h . the moment equa- 
tions are hyperbolic, which is important because otherwise the radiation model would be 
physically meaningless. The main disadvantage is that the maximum entropy closure is un- 
able to give the correct Rosseland mean in the near-equilibrium limit, and can thus not be cor- 
rec t. For example, fo r o~ v = the near-equilibrium effective absorption coefficients are given 
by iStruchtrupT l ll996l) ^ 

K t' ME - r— T / V-?) 

$™ dvv*d v n\f q) 

which is a Planck-like mean that averages k v instead of averaging its inverse. It is only seem- 
ingly surprising that the maximum ent ropy closure i s wrong even close to equilibrium. This 
closure concept must fail in general, as lKohlerl Jl948h has proven that for the linearized BTE 
the entropy production rate, rather than the entropy, is the quantity that must be optimized. Both 
approaches lead of course to the correct equilibrium distribution. But the quantity responsi- 
ble for transport is the first order deviation SI V = I v — B v , which is determined by the entropy 
production and not by the entropy. Moreover, it is obvious that Eq. j26t is explicitly inde- 
pendent of the radiation-matter interaction. Consequently, the distribution resulting from 
entropy maximization cannot depend explicitly on the spectral details of k v and o~ v , which 
must be wrong in ge neral. A crit i cal di scussion of the maximum entropy production closure 
was already given by lStruchtrupl dl998h : he has shown that only a large number of moments 
generalized to higher powers in frequency up to order i/ 4 , are able to reproduce the correct 
result in the weak nonequilibrium case. Consequently, despite of its ostensible mathematical 
advantages, we propose to reject the maximum entropy closure for the moment expansion of 
radiative heat transfer. A physically superior method based on the entropy production rate 
will be discussed in the next subsection. 



4.4 Minimum Entropy Production Rate Closure 

As mentioned, Kohler 1 1948) has proven that a minimum entropy production rate principle 
holds for the linearized BTE. The appl ication of this principle to moment expansi ons has been 
shown by lChristen & Kassubekl {2009h for the photon gas and by lChristenl |2O10h for a gas of 
independent electrons. The formal procedure is fully analogous to the maximum entropy clo- 
sure, but the functional to be minimized is in this case the total entropy production rate, which 
consist of two parts associated with the radiation field, i.e., the photon gas, and with the LTE 
matter. The latter acts as a thermal equilibrium bath. The two success factors of the applica- 
tion of this closure to radiative transfer are first that the RTE is linear not only near equilibrium 
but in the whole range of l v (or f v ) values, and secondly that the entropy expression Eq. 1 1261 

5 Convexity refers here to the mathematical entropy definition with a sign different from Eq. l26l 
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is valid also far from equilibrium (cf. lLandau & Lifshitzl te005|) '). 

In order to derive the expression for the entropy production rate, S, one can consider sepa- 
rately the two p artial (and spatial ly local) rates S rat j and S m of the radiation and the medium, 
respectively (cf. IStruchtrupl jl998h ). S rac ( is obtained from the time-derivative of Eq. l l26t , use 
of Eq. 0, and writing the result in the form 9(S rac j + V ■ J s = S rac j with 



(30) 



where n v is given by Eq. i27t . Js is the entropy current density, which is not of further interest 
in the following. The entropy production rate of the LTE matter, S ma t, can be derived from the 
fact that the matter can be considered locally as an equilibrium bath with temperature T(x). 
Energy conservation implies that W in Eq. (|8j is related to the radiation power density in Eq. 
fi"5t by W = — Pjd. The entropy production rate (associated with radiation) in the local heat 

bath is thus S ma t = W/T = —P^/T. Equation (3) implies hv /k%T — ln(l + 1/njf 1 "), and one 
obtains 

r i A +n ( e ?) \ 

Smatllv] = -* B y dvdO— In ^ 1 £(B„ - l v ) . (31) 
The total entropy production rate S = S rac j + S ma t is 

S[l v ] = dvS v = -k B f dvdd^-ln ( n ^ l + n v q) ) \ £( 
J hv \4 eq '(l + n v )J 



The closure receipt prescribes to minimize S[I V ] by varying I v subject to the constraints that 
the moments E, F, ... etc. are fixed. The solution l v of this constrained optimization problem 
depends on the values E, F, ... . The number N of moments to be taken into account is in 
principle arbitrary, but we still restrict the discussion to E and F. After introducing Lagrange 
parameters Ag and Ap, one has to solve 

S[Iv] -A E \E-^ J dvdCllJ) -Ap- (f- ^ J dvdClClIv^j =0 , (33) 

where Sj v denotes the variation with respect to I v . The solution of this minimization problem 
provides the nonequilibrium state I v . 



5. Effective Transport Coefficients 

We will now calculate the effective transport coefficients ici 6 ^, K f eff ', and the Eddington factor 
X with the help of the entropy production rate minimization closure. We assume F = (0,0, F) 
in ^-direction, use spherical coordinates (6,<p) in Q-space, such that l v is independent of the 
azimuth angle (p. For simplicity, we consider isotropic scattering with p(fl, O) = 1, although it 
is straightforward to consider general randomly oriented scatterers with the phase function p v 
being a series in terms of Legendre polynomials P n (f()- Here, we introduced the abbreviation 
]i = cos(0). Withdfi = 2nsin(8)d8 = —2ndpi, the linear operator C, acting on a function f v (}i), 
can be written as ^ 

C<p v =K v <p v (p)+<r v (tp v fa)- -J dfifyiflU , (34) 
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which has an eigenvalue k v with eigenfunction Pq(}i) and (degenerated) eigenvalues k v + a v 
for all higher order Legendre polynomials P„ (fi), n = 1,2, ... . In the following two subsections 
we focus first on limit cases that can be analytically solved, namely radiation near equilibrium 
(leading order in E — E^' and F), and the emission limit (leading order in E, while < F < E). 
In the remaining subsections the general behavior obtained from numerical solutions and a 
few mathematically relevant issues will be discussed. 

5.1 Radiation Near Equilibrium 

Radiation at thermodynamic equilibrium obeys I v = B v and F = 0. Near equilibrium, or weak 
nonequilibrium, refers to linear order in the deviation Sl v = I v — B v . Higher order corrections 
of the moments E = E^' + SE and F = SF are neglected. Because the stress tensor is an 
even function of 5 I v , x = 1 /3 remains still valid in the linear nonequilibrium region (except 
for the singular case of Auer's VEF with j = 1). We will now show that, in contrast to the 
entropy maximization closure, the entropy pr oduction minimization closu re yields the correct 
Rosseland radiation transport coefficients (cf. IChristen & Kassubekl |2009h l. 
For isotropic scattering it is sufficient to take into account the first two Legendre polynomials, 
1 and }i: 5l v = c[ ' + c v l \i, with ^-independent c v 0,1 ' that must be determined. Equations {9) 
and QS) yield 

5E V = ^j\^+c v %) = ^ (35) 
SF V = %£d?(& (36) 



and from Eq. 



^vi^') H°' )2+ ^ (w) (4 "' 2 )- 

Minimization of S v with respect to c v ° ' 1 ' under constraints SE = J dvSE v and SF = J dv5F v 
leads to 

c v = r^° E ( 3 °) 

inK v Jdvv i K v - 1 d v n { J t,> 

C (D = 3cv%n^ gp ^ (39) 

4tt(% + a v ) J dvv^iKy + a v )- l d v n v eq) 

where we made use of the relation d v n.y = n v (1 + ny e ^)h/kg T. As SI V is known to leading 
order in SE and SF, the transport coefficients can be calculated. One finds 

(eff) 2n f C{S1 V ) 4n f c[ 0) 

(eff) In f , , C(Sl v ) 47T f , , . Cy 1 ' , , 

4 ——J dvd FF £p = — / dv ( K v + "vJ^gF = \ Kv + ^/Ro ' ( 41 ) 

hence the effective absorption coefficients are given by the Rosseland averages Eqs. j2H and 
(22) . Similarly, it is shown that Ely = (E/3)<%;. This proves that the minimum entropy pro- 
duction rate closure provides the correct radiative transport coefficients near equilibrium. 
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5.2 Emission limit 

While the result of the previous subsection was expected due to the general proof by lKohlerl 
the emission limit is another analytically treatable case, which is, however, far from 
equilibrium. It is characterized by a photon density much smaller than the equilibrium den- 
sity, hence l v <C B v , i.e., E <C i.e., emission strongly predominates absorption. To leading 
order in fly, the entropy production rate becomes 

S v = -2nk B [\ dptdkhxinv) (42) 



such that constrained optimization gives 



h = ~~ VT > (43) 

with Lagrange parameters Ae and Ay. The /.(-integration in Eqs. (9) and (10) can be performed 
analytically, yielding 

E = WM ln / A £ + A F \ (44) 



where we introduced 



c 2 Af \ Aj; — Af 

T[k v ) = 47T r dvv 2 K v n v eq) . (46) 
Jo 

Up to leading order in l v , one finds by performing the integration analogous to Eqs. d40t and 

m 

4^ = <^> H and K (e ff) = TMK (47) 

/ \K V ) 

As one expects, in the emission limit the effective absorption coefficients are Planck-like, i.e., 
a direct average rather than an average of the inverse rates like Rosseland averages. The 
Eddington factor can be obtained from 1133 = %E by calculating 



2tj- rco A 
II33 = — J dv J l v , 



(48) 



which leads to 

X (v) = -^v, (49) 
A f 

where the ratio of the Lagrange parameters, and thus also the VEF, depends only on v = F/E. 
This can be seen if one divides Eq. {44) by ((45). For small v, the expansion of Eqs. {44) and 
(45) gives Ag/Af = — l/3v, in accordance with the isotropic limit. In the free streaming limit 
v — > 1 from below, it holds Af — > — Ag, which follows from ln(Z) =2 — Agln^/Ap with 
Z = (Ag + Af)/ (Ae — Af ) obtained from equalizing (44) with (45) . 

For arbitrary v the Eddington factor in the emission limit can easily be numerically calculated 
by division of Eq. (44} by Eq. {45}, and parameterizing v and x with Ap/Af. The result will be 
shown below in Fig.|4]a). It turns out that the difference to other VEFs often used in literature 
is quantitatively weak. 
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While lChristen & Kassubekl j2009h disregarded scattering, it is here included. For strong scat- 
tering a v 2> k v , Eq. ((47) implies that the effective absorption coefficient K^ eff ' of the radiation 
flux is given by a special average of cr v where K v enters in the weight function. In particular, for 
frequencies where k v vanishes, there is no elastic scattering contribution to the average in this 
limit. This can be understood by the absence of photons with this frequency in the emission 
limit. 

5.3 General Nonequilibrium Case 

The purpose of this subsection is to illustrate how the entropy production rate closure treats 
strong nonequilibrium away from the just discussed limit cases. For convenience, we intro- 
duce the dimensionless frequency £ = hv/kgT. First, we consider gray-matter (frequency 
independent k v = k) without scattering (cr v = 0). In Fig. Q]a) the quantity £ 3 «, being propor- 
tional to I v , is plotted as a function of f for F = and three values of E, namely E = E^ cc l\ 
E = E™ /2, and E = 2E^ ec l\ The first case corresponds the thermal equilibrium with l v = B v , 
while the others must have nonequilibrium populations of photon states. The results show 
that the energy unbalance is mainly due to under- and overpopulation, respectively, and only 
to a small extent due to a shift of the frequency maximum. 

Now, consider a non-gray medium, still without scattering, but with a frequency dependent 
k v as follows: K = 2%\ for f < 4, with constant K\, and k = K\ for f > 4. The important property 
is that k v is larger at low frequencies and smaller at high frequencies. The resulting distribu- 
tion function, in terms of £ 3 «, is shown in Fig.Q]b). For E = EW, the resulting distribution is 
of course still the Planck equilibrium distribution. However, for larger (smaller) energy den- 
sity the radiation density differs from the gray-matter case. In particular, the distribution is 
directly influenced by the K v ,-spectrum. This behavior is not possible if one applies the max- 
imum entropy closure in the same framework of a single-band moment approximation. A 
qualitative explanation of such behavior is as follows. Equilibration of the photon gas is only 
possible via the interaction with matter. In frequency bands where the interaction strength, 
K v , is larger (£ < 4), the nonequilibrium distribution is pulled closer to the equilibrium dis- 
tribution than for frequencies with smaller k v . This simple argument explains qualitatively 
the principal behavior associated with entropy production rate principles: the strength of the 
irreversible processes determines the distance from thermal equilibrium in the presence of a 
stationary constraint pushing a system out of equilibrium. 
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\ % 
Fig. 1. Nonequilibrium distribution (£ 3 n v oc I v ) as a function of f = hv /k%T, without scattering, 
for F = and £ = E^) (solid), E = E( e 9)/2 (dashed), and E = 2EW (dotted), a) gray matter; 
b) piecewise constant k with k^ <4 = 2k> >4 . 

Results for the effective absorption coefficients k^ S ' and Kp are shown in Fig. 12 In Fig.|2]a) 
it is shown that the effective absorption coefficient fci is equal to the Planck mean (1.26k 17 
dashed-double-dotted) in the emission limit E/EW — > 0, and equal to the Rosseland mean 
(1.6 jq, dashed-dotted) near equilibrium E = E^), and eventually goes slowly to the high fre- 
quency value K\ for large E. The effective absorption coefficient obtained from the maximum 
entropy closure is also plotted (dotted curve), and although correct for E / E ( e 1> — > 0, it is wrong 
at equilibrium E = E^\ For the present example the maximum entropy closure is strongly 
overestimating the values of K^ eff ' . 

Figure |2]b) shows K^' as a function v, for various values of E. As at constant E, increas- 
ing v corresponds to a shift of the distribution towards higher frequencies in direction of F, a 

decrease of kL must be expected, which is clearly observed in the figure. 
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Fig. 2. a) Effective absorption coefficients for £ as a function of E for F — 0, with the same 
spectrum as for Fig. Q]b). Dashed-dotted: Rosseland mean; dashed-double-dotted: Planck 
mean; solid: entropy production rate closure (correct at E — EW'); dotted: entropy closure 
(wrong at E = E^). b) Effective absorption coefficients for F as a function of v = F/E for 
different E-values (dotted: E/EW = 2; solid E/EM = 1; dashed: E/EM = 0.5; short-long 
dashed: E/E^ = 0.05). Dashed-dotted and dashed-double dotted as in a). 

In order to investigate the effect of scattering a v ^ 0, we consider the example of gray absorb- 
ing matter, i.e., constant k v = K\, having a frequency dependent scattering rate p><4 = and 
cj>4 = Kj. Scattering is only active for large frequencies. The distribution ^n v of radiation 

with E = 2E^\ with finite flux v = 0.25 for different directions }i = cos(6) = -1, -0.5, 0, 0.5, 1 
is plotted in Fig.|3]a). Since the total energy of the photon gas is twice the equilibrium energy, 
the curves are centered around about twice the equilibrium distribution. As one expects, the 
states in forward direction (pi = 1) have the highest population, while the states propagating 
against the mean flux (ji = — 1) have lowest population. This behavior occurs, of course, also 
in the absence of scattering. One observes that scattering acts to decrease the anisotropy of the 
distribution, as for £ > 4 the curves are pulled towards the state with yi « 0. Hence, also the 
effect of elastic scattering to the distribution function can be understood in the framework of 
the entropy production, namely by the tendency to push the state towards equilibrium with a 
strength related to the interaction with the LTE matter. 
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% E/E?" 
Fig. 3. a) Nonequilibrium distribution (£ 3 n v <x l v ) as a function of £ = hv/kgT, for a medium 
with constant absorption k v = K\ and piecewise constant scattering c^,, with p> <4 = 0, and 
Cg>4 = Ki. The different curves refer to different radiation directions of ]i = — 1, —0.5, 0, 0.5, 1 
(solid curves in ascending order) from photons counter-propagating to the mean drift F to 

photons in F-direction. b) Effective absorption coefficients K^ eff ' as a function of E/E^ for 

v = 0.125, 0.25, 0.5 (solid curves in ascending order); dashed-dotted: Rosseland mean, dashed: 

■ ■ c ( eff ) 

emission mean of k v . 



The effective absorption coefficient k4 is shown in Fig. |3]b) for various values of v; it is 
obvious that it must increase for increasing v and for increasing E. The Rosseland and Planck 

averages of k v + a v are given by 1.42kj and 1.40 jq, while the emission limit for Kp 6ff ' given in 
Eq. @3 is 1.20ki. 

The VEF will be discussed separately in the following subsection, because its behavior has not 
only quantitative physical, but also important qualitative mathematical consequences. 

5.4 The Variable Eddington Factor and Critical Points 

A detailed discussio n of general mathematical properties and conventional closures is given 
by lLevermord (1996). A necessary condition for a closure method is existence and uniqueness 
of the solution. It is well-known that convexity of a minimization problem is a crucial 
property in this context. One should note that convexity of the entropy production rate in 
nonequilibrium situations is often intr oduced as a presu mption for further considerations 
rather t han it is a proven property ( cf . iMartyushevI ^20061) 1. For the case without scattering, 
cr v = 0, IChristen & Kassubekl J2009|) have shown that the entropy production rate d33t is 
strictly convex. A discussion of convexity for a finite scattering rate goes beyond the purpose 
of this chapter. 

Besides uniqueness of the solution, the moment equations should be of hyperbolic type, in 
order to come up with a physically reasonable radiation model. It is an advantage of the 
entropy maximization closure that uniquene ss and hyperbolic ity are fulfilled and are related 
to the convexity properties of the entropy (cf . iLevermor el jl996k In the following, we provide 
some basics needed for understanding the problematic of hyperbolicity, its relation to the VEF 
and the occurrence of critical points. The latter is practically relevant because it affects the 
modelling of the boundary c onditions, particularly in the context of nu merical simulation s. 
More details are provided by lKorner & Tankal dl992USmit et all dl997l) , and lPons et"al] teOOOl) . 
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A list of the properties that a reasonable VEF must have (cf. IPomrariingl jl982h ) is: x(v = 
0) = 1/3, x{ v = 1) = 1/ monotonously increasing x( v )> an d the Schwarz inequality v 2 < x{ v )- 
The latter follows from the fact that x and v can be understood as averages of }i 2 and \i, re- 
spectively, with (positive) probability density l v {\i) IE. Hyperbolicity adds a further require- 
ment to the list. Equations d!3t and | |141 form a set of quasilinear first order differential equa- 
tions. For simplicity, we consider a one dimensional position spac^l with coordinate x with 
< x < L, and variables E > and F. In this case we redefine F, such that it can have either 
sign, — E < F < E. We assume flux in positive direction, F > 0, write the moment equations in 
the form 



C ' V f / V 3 E^ £ ) E3 F* J \ F J \ P F 

For spatially constant E and F, small disturbances of SE and SF must propagate with well- 
defined speed, implying real characteristic velocities. Those are given by the eigenvalues of 
the matrix that appears in the second term on the left hand side of Eq. d50t and which we 
denote by M: 

TrM / (TrM) 2 
w± = ^± ] j^fL-det(M) , (51) 

where "Tr" and "det" denote trace and determinant. Note that the w± are normalized 
to c, i.e. — 1 < W- < w + < 1 must hold. Hyperbolicity refers to real eigenvalues w± and 
to the existence of two independent eigenvectors. The condition for hyperbolicity reads 

(d F ( X E)) 2 +4d E (xE)>0. 

Provided hyperbolicity is guaranteed, the sign of the velocities is an issue relevant for the 
boundary conditions. Indeed, the boundary condition, say at x = L, can only have an effect 
on the state in the domain if at least one of the characteristic velocities is negative. It is clear 
that a disturbance near equilibrium (v = 0) propagates in ±x direction since w + = —W- due 
to mirror symmetry. Hence W- < < w+ for sufficiently small v. In this case boundary 
conditions to both boundaries x = and x = L have to be applied as in a usual boundary 
value problem. However, for finite v, reflection symmetry is broken and w + ^ —W-. It 
turns out, that for sufficiently large v, either w+ or W- can change sign. For positive F, we 
denote the value of v where W- becomes positive by v c . This is called a critical point because 
det(M) = w + W- vanishes there. Beyond the critical point, all disturbances will propagate in 
positive direction, and a boundary condition at x = L is not to be applied. This can introduce 
a problem in numerical simulations with fixed predefined boundary conditions. The rough 
physical meaning of the critical point is the cross-over from diffusion dominated to streaming 
dominated radiation. In the latter region it might be reasonable to improve the radiation 
model by involving higher order moments or partial moments, for example by decomposing 
the moment s in backward and forward propagating components E± and F± (cf. sect. 3.1 in 
iFrank! l|20u7h l 

In Fig. H] a), different VEFs are shown. All of them exhibit the above mentioned properties, 
x{v = 0) = 1/3, monotonous increase, x( v ~ * 1) — 1/ an d the Schwarz inequality v 2 < x- 
In particular, the VEFs obtained from entropy production rate minimization is shown for 
E = E^ 1 ?' for gray matter with a v = 0, as well as for the emission limit (cf. Eqs. f44t and 
(45)). Note that the latter x(v) is a function of v only and is independent of the detailed 
properties of the absorption and scattering spectra. The similarity of the differently defined 



6 Momentum space remains three dimensional. 
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VEFs, combined with the error done anyhow by the two-moment approximation, makes it 
obvious that for practical purpose the simple Kershaw VEF (j = 2) may serve as a sufficient 
approximation. 

In Fig. E]b) the characteristic velocities w± are plotted versus v for the various VEFs discussed 
above. It turns out that the VEF given by Eq. (25) has a critical point for j > 3/2 given 
by v c = l/^/2(/ — 1), and that there is a minimum v c value of 0.63 at = 3.16. The VEF 
by Kershaw and maximum entropy have v c = 1/ \/2 and v c = 2\/3/5, respectively. Also 
the VEF associated with the entropy production rate has generally a critical point, which 
depends on E. One has to expect a typical value of v c ~ 2/3. For the VEF l(25) with ; = 1 a 
critical point will not appear. In the framework of numerical simulations, this advantage can 
outweigh in certain situations the disadvantage of the erroneous anisotropy in the v — > limit. 




v v 
Fig. 4. a) Eddington Factors x versus v and b) characteristic velocities w± (see Sect. 15. 4t for 
various cases. Minimum entropy production: E = (thick solid curve) and emission limit 
E <C E^> (thin solid curve); maximum entropy (dashed); Kershaw (dotted; j = 2 in Eq. (25\), 
and Auer (dashed-dotted; j = 1 in Eq. {25)). 

6. Boundary Conditions 

In order to solve the moment equations, initial and boundary conditions are required. While 
the definition of initial conditions are usually unproblematic, the definition of boundary con- 
ditions is not straight-forward and deserves some remarks. In the sequel we will consider 
boundaries where the characteristic velocities are such that boundary conditions are needed. 
But note that the other case where boundary conditions are obsolete can also be important, 
for example in stellar physics where, beyond a certain distance from a star, freely streaming 
radiation completely escapes into the vacuum. 

The mathematically general boundary condition for the two-moment model is of the form 

aE + bn ■ F = I , (52) 

with the surface normal A, and where the coefficients a, b, and the inhomogeneity T 
must be determined in pri nciple from Eq.((5). There is a certain ambiguity to do this (cf. 
iDuderstadt & Martini {197^ 1 and thus a number of different boundary conditions exist in the 
literature (cf.lSul feOOOh ). 

There may be simple cases where one can either apply Dirichlet boundary conditions E(xvv) — 
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E w to E, where E w is the equilibrium value associated with the (local) wall temperature, 
and/ or homogeneous Neumann boundary conditions to F, (n ■ V)F = 0, at x w . This approach 
may work, if the boundaries do not significantly influence the physics in the region of interest, 
e.g., in the case where cold absorbing boundaries are far from a hot radiating object under in- 
vestigation. It can also be convenient to include in the simulation, instead of using boundary 
conditions, the solid bulk material that forms the surface, and to describe it by its K v and cr v . In 
the next section an example of this kind will be discussed. If necessary, thermal equilibrium 
boundary conditions deep inside the solid may be assumed. In this way, it is also possible to 
analytically calculate the Stefan-Boltzmann radiation law for a plane sandwich structure (hot 
solid body)-(vacuum gap)-(cold solid body), if an Eddington factor l l25t with j = 1 is used and 
the solids are thick opaque gray bodies. 

In general, however, one would like to have physically reasonable boundary conditions at 
a s urface character ized by Eq. (5). F° r engineering applications, often boundary conditions 
by IMarshakl (19471) are used. In the following, we ske t ch the principle how these boundary 
conditions can be derived (cf. iBayazitoelu & Hieenvi I ^1979^ 1. For other types, like Mark or 
modified Milne boundary conditions (cf. H3H000)). Let the coordinate x > be normal to 

the surface at x = 0, and ask for the relation between the normal flux F, E, and E^ at x = 0. 
The F-components tangential to the boundary are assumed to vanish, and diffusive reflection 
with r(x w , O, O) = r/n with r = 1 — e is considered. In terms of moments, the radiation field 
is given by 

Iv = it ( E ^o(/0 +3FvPi(F) + |(3n„,ii - E V )P 2 ( F ) + , (53) 

with Legendre polynomials Po = 1, P\ = }i, P2 = (3^ 2 — l)/2, and where Yin = The 
exact solution contains also higher order Legendre polynomials, as indicated by the dots. The 
boundary condition (5) can be written as 

I V ( F > 0) = eB v + 2r dft\fl\ I v (fi) . (54) 

By using Eq. (53), the integral can be calculated, such that the right hand side of Eq. (54) 
becomes a constant with respect to }i, while the left hand side is, according to Eq. (53), a 
function of }i defined for < y. < 1. In order to obtain the required relation between F and 
E, one has to multiply Eq. (54) with a weight function h(}i) and inte grate over u fro m to 1. 
The above mentioned ambiguity lies in the freedom of choice of h(}i). IMarshakl Jl947l) selected 
h = P\. Provided P n for n > 3 are neglected in Eq. (53), the integration leads to 



2(2- 



(3E + 15nj 



(55) 



If higher order moments are to be cons idered, additional projections h ave to be performed, in 
analogy to the procedure reported by iBayazitoglu & Higenyi I (19791) for the P-3 approxima- 
tionQFor isotropic radiation with x = 1/3, or Flu = E/3, the pref actor of E becomes unity 
and Eq. (55) reduces to the well-known P-l-Marshak boundary condition. In the transparent 
limit with x = 1/ the prefactor becomes 9/4. 

For the simple case of two parallel plane plates (e = 1) with temperatures associated with 



7 Note that ne ither the series 1531 stops after the N'th moment (even not for the P-N approximation, 
cf. ICullenl <200ll» . nor all higher order coefficients drop out after projection of Eq. {54} on P„. A general 
discussion, however, goes beyond this chapter in will be published elsewhere. 
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£ W( i and E w 2 < E w i, and separated by a vacuum gap, both moments E and F are spatially 

constant and the Stefan-Boltzmann law F = (E^ — Ei; )/4 is recovered. But note that the 

energy density E between the plates is not equal to the expected average of e[ cc,) and E^ eq) , 
which is an artifact of the two-moment approximation with VEF. 

7. A Simulation Example: Electric Arc Radiation 

The two-moment approximation will now be illustrated for the example of an electric arc. The 
extreme complexity of the full radiation hydrodynamics is obvious. Besides transonic and tur- 
bulent gas dynamics, which is likely supplemented with side effects like mass ablation and 
electrode erosion, a temperature range between room temperature and up to 3C/000 K is cov- 
ered. In this range extremely complicated absorption spectra including all kinds of transitions 
occur, and the radiation is far from equilibrium although the plasma can often be considered 
at LTE. Last but not least, the geometries are usually of complicated three-dimensional na- 
ture wit hout much symmetry, as for instance in a elec t ric circuit break e r. More details ar e 
given bvlTones & Fand <198Ch,lAubrecht & Lowkd Jl994l) , lEby et alj Jl998l) , lGodinet al.1 feOOCh , 



iDixon et al. I <2004 . and lNordborg & Iordanidisl <2008l) 
It is sufficient for our purpose to restrict the considerations to the radiation part for a given 
temperature profile, for instance of a cylindrical electric arc in a gas in front of a plate with 
a slit (see Fig. |5j. We may neglect scattering in the gas (cr v = 0) and mention that an elec- 
tric arc consists of a very hot, emitting but transparent core surrounded by a cold gas, which 
is opaque for some frequencies and transparent for others. First, one has to determine the 

effective transport coefficients fcl , x4 , and xi v )> with the above introduced entropy pro- 
duction minimization method. For simplicity, we assume now that this is done and these 
functions are given simply by constant values listed in the caption of Fig. [5] and that x(v) is 
well-approximated by Kershaw's VEF. Note that due to the low density in the hot arc core, the 
effective absorption coefficient there is smaller than in the surrounding cold gas. Therefore, 
one expects that the radiation in the arc center will exhibit stronger nonequilibrium than in 
the surrounding colder gas. 

The energy density E and the velocity vectors v = F/E obtained by a simulation with the 
commercial software ANSYS® FLUENT® are shown in Figs. |5] At the outer boundaries, 
homogeneous Neumann boundary conditions are used for all quantities. The wall defining 
the slit is modelled as a material with either a) high absorption coefficient or b) high scattering 
coefficient. The behavior of the velocity vector field clearly reflects these different boundary 
properties. The E-surface plot shows the shadowing effect of the wall when the arc radiation 
is focused through the slit. The energy densities E along the x-axis are shown in Fig. [6] a) for 
the two cases. One observes the enhanced E in the region of the slit for the scattering wall. 
The energy flux in physical units, i.e., cF , on the screen in front of the slit is shown in Fig. [6] 
b). The effect here is again what one expects: an enhanced and less focused power flux due to 
the absence of absorption in the constricting wall. 
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Fig. 5. Illustrative simulations of the moment equations with FLUENT® for a cylindrical elec- 
trical arc (radius 1 cm, temperature 10' 000 K, k£ 



eff) 
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1/m) in a gas (ambient tempera- 



ture 300 K, = K { p U> =5/m). A solid wall (a): only absorbing with x4' 

b) wall with scattering coefficient Jci = 5/m, k4 = 500/wj) with a slit in front of the arc 
focuses the radiation towards a wall. Surface plot for E (dark: large, bright: small, logarithmic 
scale); arrows for v (not F!). Only one quadrant of the symmetric arrangement is show. 
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Kt eff) = 500/m; 




8. Summary and Conclusion 

After a short general overview on radiative heat transfer, this chapter has focused on truncated 
moment expansions of the RTE for radiation modelling. One reason for a preference of a 
moment based description is the occurrence of the moments directly in the hydrodynamic 
equations for the matter, and the equivalence of the type of hyperbolic partial differential 
equations for radiation and matter, which allows to set numerical simulations on an equal 
footing. 
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The truncation of the moment expansion requires a closure prescription, which determines the 
unknown transport coefficients and provides the nonequilibrium distribution as a function 
of the moments. It was the main goal of this chapter to introduce the minimum entropy 
production rate closure, and to illustrate with the help of the two-moment approximation that 
this closure is the one to be favored due to the following properties of the result: 

• It is exact near thermodynamic equilibrium, and particularly leads to the Rosseland 
mean absorption coefficients. 

• It exhibits the required flux limiting behavior by yielding reasonable variable Eddington 
factors. 

• It gives the expected results in the emission limit, and particularly leads to the Planck 
mean absorption coefficient. 

• It can be generalized to an arbitrary number of and type moments. 

• It can be generalized to particles with arbitrary type of energy-momentum dispersion 
(e.g. massive particles) and statistics (Bosons and Fermions), as long as they are de- 
scribed by a linear BTE. In stellar physics, for instance, neutrons or even neutrinos can 
be included in the analogous way. 

The main requirement of applicability is that the particles be independent, i.e., they interact 
on the microscopic scale only with a heat bath but not among each other. On a macroscopic 
scale, long-range interaction (e.g., Coulomb interaction) via a mean field may be included on 
the hydrodynamic level of the moment equations. Independency, i.e. linearity of underlying 
Boltzmann equation, has the effect that on the level of the BTE (or RTE) nonequilibrium is 
always in the linear response regime. In this sense, all transport steady-states are near equi- 
librium eve n if f v strongly deviates from fy , and the entropy production rate optimization 
according to lKohlerl 1 19481) can be applied. 
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